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We determine the ground-state properties of a gas of interacting bosonic atoms in a one-dimensional optical 
lattice. The system is modelled by the Bose-Hubbard Hamiltonian. We show how to apply the time-evolving 
block decimation method to systems with periodic boundary conditions, and employ it as a reference to find the 
ground state of the Bose-Hubbard model. Results are compared with recently proposed approximate methods, 
such as Hartree-Fock-Bogoliubov (HFB) theories generalised for strong interactions and the variational Bijl- 
Dingle-Jastrow method. We find that all HFB methods do not bring any improvement to the Bogoliubov theory 
and therefore provide correct results only in the weakly-interacting limit, where the system is deeply in the 
superfluid regime. On the other hand, the variational Bijl-Dingle-Jastrow method is applicable for much stronger 
interactions, but is essentially limited to the superfluid regime as it reproduces the superfluid-Mott insulator 
transition only qualitatively. 



PACS number: 03.75.Hh, 03.75.Lm, 05.30.Jp 
Keywords: optical lattice, superfluid-to-Mott insulator tran- 
sition, Hartree-Fock-Bogoliubov, Bijl-Dingle-Jastrow, TEBD 



I. INTRODUCTION 

Optical lattices loaded with cold atomic gases have provided 
a fertile experimental testing ground for the study of fundamen- 
tal phenomena exhibited by quantum degenerate gases |H1|2]|3]. 
One of the greatest advantages of these systems is that mi- 
croscopic theories can be directly compared with experiments 
without the use of any fitting parameters thanks to the dilute- 
ness of atomic gases and the flexible and precise controllability 
of parameters. In particular, the Bose-Hubbard model (4] [5) 
has been successful in treating systems of cold bosonic atoms 
confined in optical lattices and explaining their various intrigu- 
ing physical properties observed in experiments, such as the 
superfluid-Mott insulator transition HHTIE), the superfluid crit- 
ical velocities |9l[l0], and the quantum depletion of conden- 
sates (TTJH2. 

The one-dimensional (ID) Bose-Hubbard model has the fol- 
lowing form 

M jj M 

Hbh = -^£(ajfl/+l+h.c.) + — £n/(n;-l) (1) 

1=1 1 1=1 

where «/ (a,) creates a boson at the lowest level localised on the 
/-th site of a ID lattice, and «; = at a; is the number operator. 
J is the hopping energy from one site to the nearest neighbor, 
and U > is the onsite repulsion. When the total number of 
particles N is incommensurate with the number of sites M, the 
system is in a superfluid state at zero temperature. When N is 
commensurate with M, the system is superfluid for small U/J 
and becomes a Mott insulator for large U/J. 

In the superfluid phase, several approximations are available 
to treat this Hamiltonian analytically. The most famous one is 
the Bogoliubov approximation [13|. However, it is limited to 
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very weak interactions, i.e. U/(nJ) <C 1, where n = N/M is 
the filling factor. It would be desirable to extend this approxi- 
mation to stronger interactions. In the context of magnetically- 
trapped ultracold atomic gases, several extensions of the Bo- 
goliubov approximation have been proposed in order to ac- 
curately reproduce the condensate fraction and the collective 
mode frequencies observed in experiments. These extensions 
are based on the Hartree-Fock-Bogoliubov (HFB) approxima- 
tion [14 1, which is closely related to the Girardeau- Arnowitt 
pair theory [38 1 . The most used one is the HFB-Shohno ap- 
proximation [16] (also known as HFB-Popov 11171 "), However, 
in this approximation, the condensate-condensate correlation is 
removed by hand to ensure a gapless excitation spectrum. Sev- 
eral authors have suggested an improved HFB approximation 
by introducing that correlation by hand so that a gapless exci- 
tation spectrum is preserved [18 19 1. More recently, Yukalov 
et al. [20 21 1 proposed a treatment of the HFB approximation 
in the context of representative statistical ensembles which pre- 
serves both the correlation and absence of a gap, and aims at 
describing the strongly-interacting regime. 

Another line of approach at zero temperature is the varia- 
tional method, which consists in minimizing the energy within 
a given subspace of the Hilbert space. The accuracy and com- 
plexity of the method depend on the choice of the subspace. 
Recently, the Bijl-Dingle-Jastrow form Il22ll23l was proposed 
as an interesting variational ansatz for bosons in a lattice 11241 . 

The purpose of this paper is to examine the "performance" 
of these theories in the case of a ID lattice system. An advan- 
tage of lattice systems is that they eliminate issues and ambigu- 
ities associated with ultraviolet divergences found in continu- 
ous systems with contact interactions. We choose a ID system 
because low dimensionality increases the effects of quantum 
fluctuations and therefore the possible differences between the 
theories. Moreover, for ID systems, we can obtain quasi-exact 
numerical solutions employing the time-evolving block deci- 
mation (TEBD) method 11251 . which can be used as a refer- 
ence. To put some perspective about the computational effort 
required by these methods, we should note that Bogoliubov- 
related calculations typically take several seconds, variational 
Bijl-Dingle-Jastrow calculations several minutes, and TEBD 
calculations several days[41 1. 

The paper is organised as follows: we first explain the TEBD 
method in Section!!]] Especially, the explanation is focused on 
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the application of the TEBD method to systems with a peri- 
odic boundary condition. The extended Bogoliubov methods 
(Bogoliubov, HFB, HFB-Shohno, improved-HFB, and HFB- 



Yukalov) are detailed in Section III The Bijl-Dingle-Jastrow 
method is presented in Section |IV| Finally, we compare all 
these methods numerically in Section [V] 
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II. THE TEBD METHOD 

Since the main purpose of this work is to examine the per- 
formance of several approximate approaches, accurate ground- 
states obtained by a quasi-exact numerical method are neces- 
sary for a reference. The TEBD method is a variant of the den- 
sity matrix renormalization group (DMRG) method 11261 [27], 
which is one of the best methods available to study ID quan- 
tum lattice systems, and provides accurate time evolutions of 
many-body wave functions and ground states via an imagi- 
nary time evolution. Although homogeneous systems with pe- 
riodic boundary conditions are better suited for the analytical 
approaches that we will use later, so far the TEBD method has 
been applied to systems with open [25 1 or infinite [28] bound- 
ary conditions. In this section, we present a detailed expla- 
nation on how to apply the TEBD algorithm to systems with 
periodic boundary conditions. 

We consider a quantum lattice system with a periodic bound- 
ary condition composed of M sites, which are labeled by index 
I, I £ {I,... ,M}. We assume that the Hamiltonian consists of 
only on-site and nearest-neighbor terms and is written as 
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respectively. In general, the number of basis configurations % 
required for convergence is of the order d L l 2 in order to express 
arbitrary states [25 1. However, since the Schmidt coefficients 
A«, decay rapidly as the index a./ increases for the ground state 
or low-lying excited states, the TEBD method with relatively 
small number of states % can be quasi-exact for these states. 

In order to compute the time evolution of a state, i.e. 
e -,i/I | v P), we split the Hamiltonian into three parts H && 
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), and H e( jg e = ^K^''^ . Subsequently, we use 

the second-order Suzuki-Trotter expansion to decompose e~' Ht 
into a product of two-site operators. When the number of lat- 
tice sites M is even, 



)> #even — Ll<m<M/2(^ 



[2m] 



MA] 



M 



where K, 



[M.M+1] _ fr[M,l] 



K. 



[l,l+l[ 



(2) 



Ki '' , reflecting the periodic boundary con- 



dition. In the case of the Bose-Hubbard model, for example, 
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correspond to the on-site interaction and the 



hopping. Spanning the Hilbert space of the whole system by a 
product of local Hilbert spaces of dimension d, a many-body 
wave function of the system can be expressed as 
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In the TEBD algorithm, coefficients c/, j 2 ,...j m are decomposed 
in a particular matrix product form as 
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The vector A<$ represents the coefficients of the Schmidt de- 
composition of I 1 ?) with respect to the bipartite splitting of the 
system into [1, . . . ,1 - 1,1] : [/+ 1,1 + 2, . . . ,M], 
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for odd M. The operators e 2 and e 2 can be ef- 
ficiently applied because they act on two neighboring T ten- 
sors [25 1. However, in order to apply e~' HeA e. eS , one needs to 
deform the tensor product Eq. (Wb to place Y^ next to rW; 
this deformation is achievable by means of the swapping tech- 
niques 1 30 1 . 

The swapping techniques allow us to exchange the po- 
sitions of two neighboring Y tensors in the tensor product 
and have been successfully applied to extend the TEBD to 
the two-legged ladder geometry [31 1 or two-component Bose 
gases ll32l . We briefly review the swapping operation intro- 
duced in Ref. [30], because it is crucial for the application of 
the TEBD to systems with periodic boundary conditions. For 
simplicity, we shall explain how to swap the positions of 
and r[' +1 l. Let us express the state in a local basis that focuses 
on sites / and / + 1 as 
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Figure 1: Schematic picture explaining how to apply e~' Hai " cS to a 
state by using the swapping techniques. Open circles represent lattice 
sites and solid lines correspond to the connections of the lattice sites 
via the nearest-neighbor interactions, (a) expresses a tensor product 
in the form of Eq. 14), In (b), rI M_1 l and D M 1 are swapped and the 
new tensors are expressed as and f . In (c), subsequently, 

ppw-2] an( j p[Af](l) swa pp e( j an( j me new tensors are expressed as 
rMP) an( j f [M-2]_ ((J) expresses a tensor product after conducting 
the swappingM — 2 times. Since rW and rM( M -2) are neighboring 
in (d), one can efficiently apply e _! ^W<5 
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In order to decompose ®a/i}a ;+1 back into the tensor product 
form as the right-hand side of Eq. (Ill, one has to reshape the 



fourth-order tensor ©a/_% +1 into a matrix ©{),,«,_,},{;,+!,«,+!} 
of the dimension d% x d%, where a/_i and a/ + i are coupled 
with ji and respectively, and conduct the singular value 
decomposition of the matrix. In the swapping process, by cou- 
pling a/_i with and 0£/ + i with ji, one reshapes the © 

tensor as ©a^+i -> ^0 {+ i,q[_i},U,(n +I }- The singular value 
decomposition of the matrix ®{; ;+ i.a ; _ 1 }.{j,.a /+1 } results in the 
swapped form of the tensor product, 



a/=l 



(12) 



The newly obtained vector corresponds to the coefficients 
of the Schmidt decomposition of I^P) with respect to the bipar- 
tition of the system into [1,. . . ,1— 1,1 + 1] : [/,/ + 2,. . . ,M], 



5 l-iJ+lh^[U+2,..M} 



«/ 



(13) 



where 



¥,[1 i-U+\]^ = £ £ ^P-i]f.['+i]y/+i 



I* 



a, 



x|4> 



)\j 



l+l) 



(14) 



I* 



[Z./+2 M]\ 

a, / 



7';=l«/+i=l 



ai+i 



(15) 



Once the swapping techniques are introduced, it is straight- 
forward to apply g^'^dge 5 t G a s j- a ]; e as indicated in Fig. |T] We 
first exchange positions of sites M — 1 and M (Fig. JTTb)) and 
next exchange positions of sites M — 2 and M (Fig. Ufc)). We 
continue conducting the swapping until the site M becomes 
next to the first site as shown in Fig.[T[d). Then, e^'^^ 5 can 
be efficiently applied. Thus, the TEBD algorithm can be ap- 
plied to systems with periodic boundary conditions by means 
of the swapping techniques. 

Using the "periodic-TEBD" method explained above, we 
calculate the ground state of the Bose-Hubbard Hamiltonian 
(1) via imaginary time evolution 



l^o) 



(16) 



where I'Prj) is an initial state. We also use the number- 
conserving version of the TEBD method, which allows sub- 
stantial speed-up of the simulations [33). In the TEBD pro- 
cedure, we choose the maximum number of atoms per site 
"max = 7 (d = 8) and retain % = 150 states in the adaptively 
selected Hilbert space. 

Having obtained the ground state, we can calculate any ob- 
servables in principle by taking the average of an operator O as 

(d) = p? t \d\v g ). 



III. EXTENDED BOGOLIUBOV METHODS 

The archetypical theory describing the weakly interacting 
Bose gas near its ground state is obtained by the Bogoli- 
ubov method lfT~3l . This method follows from the observa- 
tion that the gas should be close to a purely condensed system. 
One can therefore decompose the field operator in the form 
aj = czj + 9j, where zj is the condensate mode, and treat the 
noncondensate projection 9j as a small correction. Often, it is 
convenient to treat c as a number \/Nq, where iVo is the num- 
ber of condensed atoms, while maintaining (9j) = 0. Although 
this replacement is inexact for finite-sized systems with a fixed 
number of particles, it does not affect the results concerning ba- 
sic properties. The operator 9j is then regarded as describing 
quantum fluctuations around the classical field \/NoZj- 

The Bogoliubov approximation corresponds to treating the 
fluctuations to second order in the Hamiltonian. It results in a 
Hamiltonian which is quadratic in 9, and can be diagonalised 
in the so-called quasiparticle basis. The system can then be 
described in terms of the noncondensate density matrix ff» = 
(9j9ic) and noncondensate anomalous average = (§j§k)- 
In the ground state, for a uniform system, the condensate mode 
Zj is uniform and rim and m» are translationally invariant, so 
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Figure 2: (Colour online) Energy per particle for N = 48, M = 60 as a function of U/J. The right panel is a close-up view of the weakly- 
interacting regime. The red squares correspond to TEBD calculations, while the blue circles are obtained from the Jastrow variational method. 
The thin dashed green curve corresponds to the Bogoliubov result, the thin continuous green curve to the Popov result, while the thick dashed 
orange curve is the HFB result and the thick continuous orange curve is the improved HFB/Yukalov result. 



that we can write: 
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where n q and m q are discrete Fourier transforms of hn and m*, 
and q is the lattice wave number. The q = mode (Goldstone 
mode) is omitted, which is consistent with the orthogonality 
between the condensate and noncondensate modes. We define 
the noncondensate and anomalous densities as: 
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and it follows from average number conservation that 

no = re — re. (19) 

where «o = Nq/M is the condensate density and n = N/M is 
the total density. 



The quasiparticle diagonalisation of the Hamiltonian leads 
to the following expressions for h q , m q in the ground state: 
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One can also show that the excitation energies as a function of 
the excitation momentum are given by the quantity 



O) 2 - A 2 . 



(25) 



Interactions generally tend to decrease the condensate frac- 
tion, and as a result the Bogoliubov approximation breaks 
down for sufficiently strong interactions. Especially it is unable 
to predict the Mott-insulator transition in commensurate lat- 
tice systems [34|. Even in incommensurate systems which re- 
main superfluid, it is valid only in the weakly-interacting limit, 
which can be understood as follows. In the case of ID Bose 
gases in continuum (with no lattice), the validity of the Bogoli- 
ubov approximation is judged by the dimensionless parameter 
7= / 2 /<= 2 which expresses the competition between the aver- 



age interparticle spacing / and the healing length | 
The healing length corresponds to the de Broglie wavelength 
associated with the ID interaction energy. When y< 1, the 
wave functions of the particles are well overlapped with each 
other and the gas can be approximately described by a macro- 
scopic wave function (the condensate mode) with long-range 
phase coherence. Therefore, the Bogoliubov approximation is 
valid in this regime. As y increases, the overlap of the wave 
functions becomes smaller and the physical quantities, such as 



the ground state energy (24i and the excitation spectrum (25 i 
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Figure 3: Condensate fraction for N = 48, M = 60 as a function of U /J. The right panel is a close-up view of the weakly-interacting regime. 
The same conventions as those of Fig. 2 are used. 
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Figure 4: Excitation spectrum for N = 48, M = 60. The same conventions as those of Fig. 2 are used. The left panel corresponds to U/J = 1 and 
the right panel corresponds to U/J = 8. Note that the apparent dispersion of the blue circles are residual fluctuations due to the stochastic nature 
of the variational Monte-Carlo method. 



deviate from those of the exact calculations based on the Bethe 
ansatz [37]. In the limit of y^> 1, long-range coherence disap- 
pears and the particles behave as impenetrable objects, which 
is referred to as the Tonks-Girardeau gas l38l . In the case of the 
Bose-Hubbard model, using | =h/(m*c) B5l and / = d/n, the 
dimensionless parameter is expressed as 7= U/(2nJ), where 
m* = h 2 /(2Jd 2 ) is the effective mass and c = (2nJU) l l 2 d/Ti is 
the Bogoliubov sound speed. Thus, the Bogoliubov approxi- 
mation is valid only when U <C nJ. 

One way to improve the Bogoliubov approximation is to 
include the higher-order perturbation terms, and make a de- 
coupling approximation which reduces the Hamiltonian to a 
quadratic form again. This is the case when one assumes that 
the noncondensate particles are uncorrelated (using Wick's de- 
coupling scheme for the operator 6). One then obtains the HFB 
approximation. The structure of the equations is the same as 
above, but now 
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where Eb os refers to the expression on the right-hand side of 
Eq. (24 1 - note that it is not the Bogoliubov energy since (O q 
and Anave changed. The solution is now implicit, and Eqs. 



17|18|19|20|21|26|27[ > have to be solved self-consistently. 



While it retains the elegance of the Bogoliubov approxima- 
tion, and is supposedly applicable to more strongly-interacting 
cases, the HFB approximation suffers from a notorious prob- 
lem related to the existence of a gap in the excitation spectrum 
flTl . which is forbidden in the superfluid regime according to 
the Hugenholtz-Pines theorem [39 1. This gap can be expressed 
as \in\ q ^oE q = 2U y/—mno, and clearly comes from the anoma- 
lous average m which is supposed to be strictly negative, ac- 
cording to Eq. ( pT| >. One ad hoc solution to this problem, first 
used by Shohno |fl6ll . is to set m to zero by hand in the HFB 
equations (|17|19|20|26|27|2"8 1 - note that this is not consistent 



with Eq. (21 1 and leads to the violation of conservation laws. 
This solution has been popularised as the HFB-Popov approx- 
imation in the recent literature on magnetically trapped Bose- 
Einstein condensates ifm . 

However, the importance of the anomalous average near the 
ground state has been stressed by some authors IT8l [T9l l20l 
who devised new variations of the HFB approximation in order 
to include that average. 

In the "improved HFB" approach fl8l [T9l it is argued that 
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Figure 5: Energy for N = 60, M = 60 as a function of U/J. The right panel is a close-up view of the weakly-interacting regime. The same 
conventions as those of Fig. 2 are used. 
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Figure 6: Condensate fraction for N = 60, M = 60 as a function of U /J. The right panel is a close-up view of the weakly-interacting regime. 
The same conventions as those of Fig. 2 are used. 



the decoupling of noncondensate particles is too strong an 
assumption, resulting in noncondensate particles interacting 
through a bare r-matrix. On the other hand, the condensate 
particles interact through a full r-matrix, renormalised by the 
presence of m, which is a measure of condensate-condensate 
correlations. If one assumes that the noncondensate correla- 
tions neglected in the HFB approximation amount to some- 
thing similar to m, it seems reasonable that the bare r-matrix 
for noncondensate particles should be upgraded to the same 
renormalised r-matrix for condensate particles. This can be 
done by hand in the HFB equations by changing the sign in 



front of m in Eq. (26 1. This way, the gap disappears. 

In the representative statistical ensemble approach ll20ll2TI . 
it is argued that the number of condensate particles No deserves 
a special treatment because it appears as a new macroscopic 
quantity in the system. As a result, a new chemical potential 
Ho associated to No should be introduced. It is distinct from the 
chemical potential p. associated to the number of nonconden- 



sate particles N = N — No- This distinction changes Eq. (26 1 
to 

CO q = {n Q - m)U + (jUo - A ) + 4/ sin( || f . (29) 

If one chooses jio — A = 2mU, the HFB equations become gap- 
less. Note that this results in the same modification as in the 
"improved HFB" approach. Although the two approaches are 
conceptually different and do lead to different results in more 
general cases (such as nonlocal interactions), they yield the 
same equations in the case studied in this paper. 

In both cases, the energy is obtained by formally changing 
the sign of m: 



Ehfby = E Bu g -Unm+^-U {n 2 - m z ) 



(30) 



IV. THE VARIATIONAL BIJL-DINGLE-JASTROW 
METHOD 

We now consider a variational method. It consists in mini- 
mizing the energy of the system within a subspace of the full 



Hilbert space. The minimization is numerically possible if that 
subspace is not too large, and leads to a good approximation of 
the true ground state if it is general enough to allow the impor- 
tant features of that state. As recently proposed in Ref. [24], we 
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Figure 7: Excitation spectrum for N = 60, M = 60. The same conventions as those of Fig. 2 are used. The left panel corresponds to U/J = 2 and 
the right panel corresponds to U/J = 8. Note that the apparent dispersion of the blue circles are residual fluctuations due to the stochastic nature 
of the variational Monte-Carlo method. 



choose the subspace of wave functions having the Bijl-Dingle- 
Jastrow form [22 , 23 1: 



x V{x\...Xn) 



■Xj) 



(31) 



This is a simple form obtained by a product of some correla- 
tion function / for all pairs of bosons. The Bijl-Dingle-Jastrow 
form Eq. 



(31 



is more general than the structure of the Bo- 
goliubov ground state (which can be reproduced from Eq. pT) 
in the limit of weak interaction). Thus it can be applied suc- 
cessfully to more strongly-interacting regimes. However, this 
method gives only the ground-state properties and some exci- 
tation properties, unlike the previous methods which are also 
applicable at finite temperature. Another drawback is that the 
minimum cannot be found analytically. 

To find the minimum numerically, we use the methods de- 
scribed in Ref. (24). Namely, minimization is performed by 
means of the power method, and at each step a Monte-Carlo 
algorithm computes all the average quantities relevant to the 
minimization procedure. Once the optimal / is found, any 
observable can be obtained by taking the average of the cor- 
responding operator in the quantum state x ¥. This average is 
calculated using the Monte-Carlo algorithm. 



V. RESULTS 

We apply the methods described in the previous sections to 
solve the Bose-Hubbard model in two cases: an incommen- 
surate case N = 48, M = 60 (for which the lattice is filled to 
80%) and a commensurate case N = 60, M = 60 (for which the 
lattice is filled to 100% ). The Superfluid-Mott insulator tran- 
sition occurs only in the commensurate case. In both cases, 
we look at characteristic properties of the ground state such 
as its energy per particle E and its condensate fraction no, as 
well as the excitation spectrum E q . In the extended Bogoli- 
ubov methods, these quantities appear naturally as variables of 
the problem - see Eqs. (gSJ, ([19) and ((25]). In the TEBD and 
Bijl-Dingle-Jastrow methods, the energy is obtained by calcu- 
lating the average (H), and the condensate fraction is given 
by Nq/N, where the number of the condensate bosons A^o is 



obtained from the largest eigenvalue of the one-body density 
matrix (djd/) l29l . A very close upper bound of the excitation 
spectrum is obtained from the /-sum rule 11401 . which states 
that 

^ -2sin(^) 2 
E q <K 

where both the kinetic energy K = {—JY!u=\ (djd,+i + dt +1 d,-)) 
and the structure factor S q = jjY^u{{^fij) ~ (ni}(fij})e' q ^~^ 
can be calculated by quantum averages. 

Results for the incommensurate case are shown in Figs. [2] [3] 
and |4] One can see that there is excellent agreement between 
the Bijl-Dingle-Jastrow calculations and the TEBD results for 
any interaction strength. On the other hand, all the Bogoliubov 
methods agree with these results only in the weakly-interacting 
limit U/J < 0.2. Curiously enough, the energy in the Bo- 
goliubov approximation shows very good agreement with the 
TEBD and Bijl-Dingle-Jastrow energy, even for strong interac- 
tions U / J ^ 0.2. While this might provide some analytical in- 
sights to devise a more elaborate theory, we believe that this is 
most likely a mere coincidence, since the Bogoliubov conden- 
sate fraction is correct only for U/J < 0.2. The same applies 
for the Popov approximation which reproduces the correct con- 
densate fraction up to U/J < 2 but fails to predict the average 
energy for U /J> 0.2. 

Results for the commensurate case are shown in Figs. [5] 
[6] and [7] One can see that a quantum phase transition for 
the superfluid phase to the Mott-insulator phase occurs around 
U/J = 4. Since we are dealing with a finite-size system, the 
transition is smoothed. All the extended Bogoliubov methods 
fail to reproduce that transition, since their range of validity is 
limited to U/J < 0.2, but again, the Bogoliubov approximation 
shows a surprisingly excellent prediction of the energy over the 
whole superfluid regime. As it was found in Ref. [24], the Bijl- 
Dingle-Jastrow method does show a quantum phase transition 
and accurately predicts the energy and condensate fraction in 
the strongly-interacting limit of the Mott-insulator regime. The 
excitation spectrum in Fig. [7] also shows that a gap appears in 
the Bijl-Dingle-Jastrow method, while all the extended Bogoli- 
ubov spectra remain gapless (by construction). 

While these results seem to indicate that the Bijl-Dingle- 
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Figure 8: Gap in the excitation spectrum (E q for q = 0) as a function 
of U/J, for N = 60, M = 60. The red squares represent the TEBD 
results, and the blue circles correspond to the variational Bijl-Dingle- 
Jastrow method. 

Jastrow ansatz is able to reproduce the superfluid-to-Mott in- 
sulator transition, as claimed in [24|, it should be noted that 
it does only in a superficial way. Firstly, the transition itself 
is only approximately reproduced by the Bijl-Dingle-Jastrow 
method, which predicts a somewhat surprisingly sharper tran- 
sition than the actual one and at a different value of U /J around 
5. This fact is expected since the system is known to be 
strongly correlated near the transition, in a way which is hardly 
recovered by any simple ansatz. Secondly, although the Bijl- 
Dingle-Jastrow ansatz tends towards the exact noninteracting 
ground state when U/J — > 0, it does not tend towards the exact 
Mott insulator state when U /J Indeed, while it can repro- 



duce accurately the energy and condensate fraction, it fails to 
reproduce the excitation gap, which even departs further from 
the real one as U/J — > °°, as shown in Fig. [8] From this we con- 
clude that the Bijl-Dingle-Jastrow ansatz provides a partial de- 
scription a Mott-insulator, the usefulness of which is restricted 
to the calculation of quantities deriving from low-order corre- 
lations. 



VI. CONCLUSION 

We have performed an extensive comparison of several the- 
oretical treatments of the Bose-Hubbard Hamiltonian ground 
state, describing bosonic atoms interacting in an optical lattice. 
The TEBD method was generalised to periodic boundary con- 
ditions in order to solve the problem. It was used as a quasi- 
exact reference to compare with the results of extended Bogoli- 
ubov methods and the variational Bijl-Dingle-Jastrow ansatz. 

We showed that all the methods refining the original Bogoli- 
ubov approximation do not in fact bring any significant im- 
provement. The validity of all these methods is restricted to 
the weakly-interacting regime U/J < 0.2. On the other hand, 
the Bijl-Dingle-Jastrow ansatz proves to be a very accurate 
approximation of the ground state in the superfluid phase for 
any strength of the interaction, showing that particles in such 
a phase are essentially correlated by pair. However, it gives 
only a partial account of the ground state in the Mott insula- 
tor phase, showing that more-than-two-particle correlations are 
needed in that phase. While the structure of the Mott insulator 
in the strongly-interacting regime is known, devising a theory 
which can account for both the superfluid and Mott insulator 
phases accurately is still a challenging problem. 
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